Parsing mPFC: A Joint Meta-Analytic and Graph Theoretic Approach

1. Department of Psychological and Brain Sciences, Boston University, Boston, USA

Corresponding author:

Abstract

Valuation effects are consistently observed in medial prefrontal and posterior cingulate cortex (mPFC and PCC). The spatial extent of these effects is mostly indistinguishable from the default mode network (DMN) in existing meta-analyses. However, little is known about how valuation effects fit within the broader functional architecture of mPFC and PCC, or whether that architecture is consistent or idiosyncratic across individuals. Here we complement a meta-analysis with fMRI-based graph theoretic approaches to subdivide mPFC and PCC at the single-subject level. Our results suggest the functional topography of mPFC has substantial variability across individuals. This highlights the potential usefulness of estimating brain effects at the individual level in this region, and points to limitations of aggregative methods such as coordinate-based meta-analysis in determining whether valuation and DMN effects emerge from common or separable brain systems. Our approach shows promise in addressing this issue through future manipulations of valuation.

Keywords: brain networks; DMN; individual differences

Introduction

Studies of decision making have consistently associated activity in ventromedial prefrontal cortex (vmPFC) with the subjective value (SV) assigned to decision outcomes (Bartra, McGuire, and Kable 2013; Clithero and Rangel 2014; Hiser and Koenigs 2018; Kable and Glimcher 2007; Levy et al. 2011). However, the distributed cortical valuation system also includes posterior cingulate (PCC) and more extended mPFC regions, which collectively show significant spatial overlap with the default mode network (DMN) along the medial wall (M. D. Fox et al. 2005; Laird et al. 2009; B. T. Yeo et al. 2011). This shared neural coverage has prompted the question of whether (and to what extent) these two systems can be dissociated.

Previous meta-analytic work has provided important insights on the psychological phenomena attributed to mPFC subregions. Independent coordinate-based meta-analyses have found mirroring spatial coverage in mPFC for valuation (Bartra, McGuire, and Kable 2013) and DMN (Laird et al. 2009). Whole brain analyses have shown that while more ventral mPFC regions are increasingly associated with reward-related elements of decision making, vmPFC strongly coactivates with DMN regions across the literature, with no clear topographical distinction between these phenomena (Vega et al. 2016). More focused meta-analytic work aimed at disentangling DMN from subjective value has determined that these systems are inseparable within mPFC (Acikalin, Gorgolewski, and Poldrack 2017). These findings highlight the multifaceted embedding of psychological constructs in subregions of mPFC, and suggest that DMN and valuation might indeed be subserved by the same system. However, recent activation-based meta-analytic work has demonstrated considerable domain specificity in this area (Kragel et al. 2018). While informative, a key limitation in distinguishing constructs through this approach is that meta-analyses often rely on information derived from group-averaged data. Averaging has traditionally been beneficial in identifying tendencies in brain function and organization when faced with short scanning sessions per subject. However, it has been shown that mPFC is subject to considerable idiosyncratic cortical folding (Zilles, Palomero-Gallagher, and Amunts 2013) and inter-subject functional variability (Mueller et al. 2013) compared to the rest of the brain, thus making group-averaging studies difficult to interpret.

Seeking to remedy these shortcomings, recent work has prescribed relevance to the analysis of single subjects in fMRI. Individual analyses of fMRI data have identified idiosyncratic, reliable, and valid functional organization that would otherwise be blurred in aggregative estimations (T. O. Laumann et al. 2015; Gordon et al. 2017; Gratton et al. 2018). Furthermore, subject-specific network arrangements have been found to predict behavioral characteristics (Kong et al. 2018). In regards to DMN, the trend of individualized analyses has led to finer idiosyncratic subdivisions of the DMN through careful selection of seed-based functional connectivity (Braga and Buckner 2017). It is thus possible that the indisitinguishable overlap of DMN and valuation effects can be attributed to a lack of spatial resolution that is better understood at the individual level. With this in mind, an important first step in disentangling these phenomena would be to determine the degree of topological heterogeneity of the DMN present within mPFC.

A persistent issue in analyzing mPFC in individuals is the various factors inducing signal dropout in its more ventral segments (Logothetis 2008). As such, there is value in examining all possible sources of covariation that compose DMN, instead of estimating topologies one seed at a time. Connectome-based analyses of resting state functional connectivity (rsFC) have been fruitful in characterizing individualized functional topologies that match task-induced activity (Gordon et al. 2017; T. O. Laumann et al. 2015; Tobyne et al. 2017). A popular approach to analyze these connectomes is to rely on graph theoretic methods, which provide an organic avenue to understand brain network dynamics (D. S. Bassett, Zurn, and Gold 2018). From these, community detection algorithms have been successfully utilized to section brain networks into cohesive substructures (Garcia et al. 2017). Such communities represent clusters of network nodes that are more connected with each other than with the rest of the network (Fortunato and Hric 2016). Among the algorithms used to estimate these communities, modularity has been widely effective in subsectioning brain networks into multiple groups (Garcia et al. 2017). However, in attempting to refine the DMN topology at the individual level, dividing these regions into what belongs and does not belong to DMN makes dissection through spectral partitioning (SP) a more valuable option. In short, spectral partitioning performs eigenvalue decompositon on an adjusted form of a correlation matrix (Laplacian matrix). The resulting eigenvector corresponding to the second-to-last eigenvalue (i.e. the Fiedler vector) provides values to divide the network into a positive and a negative community (Fiedler 1975; Belkin and Niyogi 2003). SP has recently been used to characterize the posterior-anterior functional gradient of the insula using rsfMRI data (Tian and Zalesky 2018). This method can thus capture a refined group of brain regions that more cohesively share activation patterns during rest.

In this study, we aim to subsection regions commonly attributed to both DMN and subjective value into subject-specific DMN and non-DMN partitions, quantifying the degree of topographical heterogeneity present in this network within and across individuals. This description will allow researchers to generate more precise topographic targets for future studies of decision making. We do this by capitalizing on the respective strengths of meta-analytic and subject-specific analyses of brain data. First, we define regions of interest (ROI) by identifying overlapping activation throughout the literature projected to an established brain atlas (M. F. Glasser et al. 2016). We then produce a rsFC networks of all the surface vertices within these ROIs for each individual resting state dataset from the Human Connectome Project (HCP; citation), and estimate the extent of DMN coverage through spectral partitioning. The resulting patterns are then evaluated for consistency across methods, throughout each individual’s time series, and across subjects. Finally, we show some interesting properties of SP related to the stability of brain networks.

Methods

Meta-analysis

We used data from metan-alyses that gathered peak brain coordinates of activity from studies of valuation (Bartra, McGuire, and Kable 2013) (n = 27 studies) and DMN (Laird et al. 2009) (n = 77). These represent the surviving areas post-statistical thresholding from each study. For each peak volumetric coordinate, 10 mm sphere masks were generated and then projected to a standard cortical mesh (fs_average, 32,000 vertices) using FreeSurfer. The sphere sections that overlapped with the cortex (should I add the details of the overlap?) were mapped to an atlas of the human cortical surface (360 regions)(M. F. Glasser et al. 2016). This produced a list of standardized parcels that were reported on each study.

With these lists, we identified the areas that were more consistently shared between DMN and valuation. For each literature, we summed the number of times each parcel was reported across studies, and performed a one-way chi-squared test of proportions on each parcel with a null probability equivalent to the bootstrapped (5000 iterations) mean of each literature’s report counts. In other words, we checked which regions were significantly (p < 0.05) more reported in their respective literatures than their counterparts. Regions of interest were then determined based on which of the remaining regions were common across literatures.

Next, we identified areas that were solely related to either literature. We limited the search space to the parcels that were contained within the DMN parcellation defined by Yeo et al. (2011). We permuted the label of every study (i.e. DMN or valuation) 5000 times while preseving the total number of studies observed in each domain, and on each iteration stored the maximum statistic from a area-wise chi-squared test of proportions on the number of permuted studies that had reported each area. This gave us a null distribution of chi-squared values that was used to determine which ROIs were more significantly prevalent on each literature.

fMRI Data

In order to quantify the intrinsic connectivity of these ROIs in the brain, we acquired resting-state fMRI data from the Human Connectome Project (HCP) (Van Essen et al. 2012) Q6 release (n = 100, randomly sampled from the total pool of 469 available subjects). Each subject’s data was acquired over two days at Washington University in St. Louis on a Siemens CONNECTOM Skyra MRI scanner (Siemens, Erlangen, Germany). Four resting state sessions (TR = 0.720 s, TE = 33.1 ms, FA = 52°, multiband factor = 8, 72 slices, 2 mm isotropic voxels) were each comprised of 1200 TRs for a total scan time of 14 min 33 s, with left-right and left-right phase encoding implemented on each day. We used this scanning structure to estimate each subject’s topology based on session, day, and overall data (4800 TRs).

Beyond the in-house minimally preprocessed pipeline from HCP (M. F. Glasser et al. 2013), which includes considerable motion correction, temporal denoising, highpass filtering (0.0005 Hz threshold), and MNI152-based normalization, all scans went through a number of additional refinements (Tobyne et al. 2017). These included band-pass filtering (allowed frequencies ranged from 0.009 and 0.08 Hz), as well as mean greyordinate signal regression (Burgess et al. 2016). Only subjects with both left-right and right-left phase encoding were included (i.e. subjects with four rsfMRI sessions). In addition, only datasets with either low motion levels (under 1.5mm) or under 0.5 mm mean framewise displacement (FD) were used. Volumes that displayed an FD of over 0.5 mm were considered as spikes and removed from the subject’s dataset, and subjects whose spike counts were above 5% of their total data were excluded from analyses. Finally, data acquired on the same day (i.e. left-right and right-left phase encoding session) were temporally demeaned. Each subject’s brain was comprised of 32k vertices per hemisphere. We retained only the cortical surfaces, which resulted in 59,412 total surface vertices per subject.

Within-ROI Network Definition

To establish each subject’s network, we selected all the vertices contained within the ROIs (n = 5,081 per subject) and computed the Pearson correlation of the time series for every pair of vertices. This produced a weighted network for each subject, where surface vertices were the nodes, and edges the correlations among all of them. Next, we applied Fisher’s r to z transformation, and took the exponential of these values so that all weights were positive while maintaining the shape of the original correlation distribution, and retained the edges whose adjusted p-values remained significant after correction for multiple tests (FDR < 0.05). As part of the evaluation step, this procedure (and the community detection outlined below) was applied at four levels: 1) each session separately (1200 TRs each); 2) the concatenated timeseries from each pair of daily sessions (2400 TRs total); 3) each subject’s full dataset (4800 TRs total); and 4) on each step of a sliding window analysis (see Evaluation for more details).

Community Detection

Communities (i.e. clusters) were identified by means of spectral partitioning (Fiedler 1975; Belkin and Niyogi 2003). First, each adjusted network was represented as an \(n\) x \(n\) correlation matrix (where \(n\) = number of vertices within the ROIs, 5,319). This matrix was then transformed into its normalized Laplacian form

\[\begin{aligned} L = I - D^{-\frac{1}{2}}WD^{-\frac{1}{2}} \end{aligned}\]

Where \(I\) is the identity matrix, \(D\) is a diagonal matrix containing the strength of each vertex (i.e. the sum of its correlations with every other vertex), and \(W\) is the correlation matrix. We then computed the eigenvalues and eigenvectors of the normalized Laplacian, and used the eigenvector associated with the second-to-lowest eigenvalue to divide the network into two. This vector (from now on called the ‘Fiedler Vector’) provides a set of positive and negative values to binarize the network with, and guarantees that the resulting communities are connected (Fiedler 1975). Given that this data-driven method is agnostic to the prevalent function of each community (i.e. DMN or not DMN), we defined the DMN community as the one containing the one DMN-specific area from the meta-analysis (ensuring that the labeling of the DMN community was consistent across subjects). Importantly, given the high density of these networks (see results), spectral partitioning was unlikely to face the issues associated with its use in sparse networks (Fortunato and Hric 2016).

In order to evaluate the validity of the resulting partitions, we also estimated these communities using the more traditional approach of modularity maximization using the algorithm from Clauset et al. (2004). This method maximizes a quality function that compares the strength of the connection between any two vertices against a null model of their probability of being randomly connected to any other vertex in the network. The method heuristically iterates through many possible combination of vertices, and selects the clustering that maximizes the network’s modularity. Thus, unlike SP, modularity can fractionate a network into many functional groups. If the partitions from the bounded and unbounded methods show high level of agreement, we can gain confidence in that we are indeed fractioning the ROIs into those belonging to DMN and not.

Partition Evaluation

In order to evaluate the stability and topographical heterogeneity of the partitions within and across subjects, we chose the adjusted rand index (ARI) as the reference metric of agreement between two clusterings (Hubert and Arabie 1985). ARI estimates the rate at which every pair of nodes in partition A were jointly assigned to the same/different parcels on partition B, over all possible changes (e.g. clustered together on A, but separately on B). Because of this, the adjusted rand index determines the percentage of agreement between any two partitions while being agnostic to the labeling scheme. The resulting ratio is compared to a baseline given by the expectation of a random node assignment across partitions. This means that even though clustering deviations are heavily penalized, positive ARI values compare favorably against chance clustering (and can take negative values if the index falls below chance).

We performed a number of comparisons among partitions. First, we computed the degree of agreement between SP and modularity per subject to validate the resulting functional patterns. Modularity and SP show overfitting and underfitting tendencies, respectively, in their community detection performance in a diverse number of network types (Ghasemian, Hosseinmardi, and Clauset 2018). Their alignment would thus give us confidence in the validity of the resulting topologies. Next, we compared the functional organization derived from the whole time series among all subjects, and report the mean RI for the group. We then estimate the level of agreement between scans taken across days (i.e. pairs of daily left-right and right-left phase encoding sessions). If the functional organization estimated by SP is indeed subject-specific, we should see higher agreement for the same subject across days than either day compared to every other subject. We then extend this idea by computing the overall agreement across all sessions. Similarly, we expect that sessions within subject will show higher agreement than between subjects. Finally, we perform a sliding window analysis (20 min windows, 1 min increments, mean number of windows per subject = 37) 8) comparing each window’s functional pattern against the clustering derived from the whole subject’s data. This approach should give us a sense of the dynamics of the functional organization of mPFC during rest. We complement this analysis by calculating the proportion of times each surface vertex was affiliated with the DMN community across time, as well as computing the flexibility of each node (i.e. the number of times each vertex changed affiliation across time). Given that the functional heterogeneity of brain areas is unevenly distributed across the brain (Mueller et al. 2013), we also performed these comparisons for PCC and mPFC separately. Given what we know about these general regions, we expected PCC to be well aligned across participants, and stable within subjects. On the other hand, we hypothesized low agreement of mPFC partitions across participants, but with considerable stability across days and sessions within subject. We compared the agreement within subjects against across subjects by computing the ratio of the mean RI within and between subjects. Ratios close to one denote good alignment across subjects, while ratios significantly higher than one suggest that partitions per subject were more stable than what you would expect against every other subject.

Given the descriptive relevance of Fiedler eigenvector values in previous fMRI research (Tian and Zalesky 2018), we tested the association of their absolute magnitude with the proportion of times a vertex was part of the DMN community during the sliding window analysis. To do this, we performed a mixed effects logistic regression model in which each subject’s eigenmap values were used as random predictors of the proportion of times they participated in the DMN community.

Estimation of Noise Effects

BOLD signal recording through fMRI in mPFC in general, and vmPFC in particular, suffers from technical issues that promote signal dropout (Logothetis 2008). It was thus important to quantify the potential influence of noise in our findings. We examined this through two complementary steps. First, we computed the temporal signal-to-noise ratio (tSNR) of all the vertices in the ROIs from HCP minimally preprocessed data (i.e. prior to the complementary denoising described in the “fMRI Data” section). We compared the mean tSNR for each region across subjects through a permutation analysis (all permutation analyses involved 5000 iterations). In addition, we compared the mean tSNR between mPFC and PCC through a permutation analysis for each separate subject, and report the percentage of subjects who had a significantly lower mean tSNR in mPFC after FDR correction for multiple testing (???). Next, we estimated which vertices within mPFC had a stronger relationship to the overall dynamics of DMN. To do this, we computed the mean time series of each subject’s estimated DMN map, which was then predicted by the time series each mPFC vertex in a mass univariate analysis. The resulting \(R^2\) values from each model were extracted to determine the total variance of the mean DMN time series explained by each individual vertex. If high noise in more ventral mPFC areas is affecting our results, we should see a gradient of \(R^2\) values that decreases ventrally across subjects. We then spatially correlated (Pearson) the \(R^2\) and tSNR values to determine whether low tSNR negatively impacted the individual contribution among vertices. (PERONSAL NOTE: THIS IS NOW DONE BY CONCATENATING ALL SUBJECTS. IT WOULD PERHAPS BE BETTER TO DO A MIXED EFFECTS MODEL WITH SUBJECTS AS RANDOM EFFECTS, AND TSNR AS A FIXED ONE)

Seed-based FC vs Community Detection

A potential downfall of the present approach is that it can be computationally demanding, while more direct seed-based functional connectivity methods could provide sufficient approximations of an individual’s DMN orgnization. We examined this issue by estimating mPFC functional maps from its vertex-wise functional coupling with the mean PCC activation. For each subject’s concatenated day 1 and day 2 scans, we computed the mean time series of all vertices contained within the PCC region of a standard DMN arrangement (B. T. Yeo et al. 2011), and correlated it (Pearons) with the time series of every mPFC vertex from our ROIs. The resulting day 1 and day 2 mPFC maps were then correlated within and between subjects in three ways: 1) between correlation-based maps; 2) between community-based maps; and 3) between correlation- and community-based maps. These spatial correlations are aimed to determine the day-to-day stability of each approach, as well as their level of agreement. These methods should produce somewhat similar arrangements, but the one displaying the highest within-subject agreement across days should be preferred (for example, see (Kong et al. 2018; Mueller et al. 2013)). We thus compared the mean within-subject inter-day correlation coefficients of each method through a paired permutation analysis (5000 iterations).

Anatomical Correspondence

Finally, we examined whether each idiosyncrasies in functional topology corresponded to outstanding anatomical features. To do this, we correlated each subject’s Fiedler vector values with their respective cortical thickness, cortical folding, and myelin density (quantified as the ratio of T1 to T2 magnitudes at every given vertex). These simple correlations are aimed to provide a coarse sense of this possible relationship, which can be studied in depth in the future.

Results

Meta-analysis

We performed a coordinate-based meta-analysis to limit the network search space to a set of overlapping and unique standard regions across DMN and SV literatures (see Methods). Volumetric coordinates from a set of DMN and SV studies were projected onto a cortical surface, and mapped to a multimodal cortical parcellation (M. F. Glasser et al. 2016) to produce a list of brain areas reported per study. Using each literature’s mean number of reports per area as the null probability (DMN = 12, 95% CI [11.53, 12.91]; SV = 5, 95% CI [4.56, 5.18]), one-way chi-squared tests of proportions were performed to identify regions that were significantly present in their respective literatures. The resulting lists were used to identify the set of common brain areas. Unique areas were then determined by permuting the affiliation of each study (i.e. DMN or SV) to create a null distribution of chi-squared statistics, which as used to identify regions more often reported in either literature. Importantly, we chose to examine these areas bilaterally, even if their significance was present in only one hemisphere.

Figure 1 shows the resulting ROIs from these analyses. DMN included 50 significant regions, while SV was comprised of 27 (see Supplemental Table 1). As expected, overlapping areas included 15 bilateral regions in PCC (31a, 31pv, RSC, d23ab) and mPFC (9m, 10r, 10v, OFC, a24, d32, p24, p32, s32), as well as a region in IPL (PF) and TPJ (PGi) (allp< 0.05). A region in visual cortex (V2) was also significant, but was removed due to its know unrelated function. Area 7m was the only one to show unique affiliation, which was significant for DMN (\(\chi^2\) = 9.12, p < 0.001).

Figure 1. Template atlas from Glasser et al. (2016) (top), as well as the common and unique brain regions associated with DMN and subjective value (bottom). Overlapping ROIs (in blue) were defined as the common brain regions that were significantly present in their respective literatures, while a single DMN-specific region (7m, in red) was identified from a permutation analysis comparing these literatures.

Partition Comparison: Full Datasets

With these 16 ROIs (15 overlapping, 1 DMN-specific), we estimated each individual’s DMN coverage using the whole time series (4800 total TRs from four scanning sessions over two days). Each subject’s network (n = 100) was constructed by performing pairwise correlations among the time series of all the surface vertices contained within the ROIs (5081 total vertices). After correcting for multiple testing, we performed spectral partitioning and modularity on these networks to identify cohesive functional communities denoting DMN and non-DMN areas. On average, networks had high densities (i.e. highly connected) regardless of the amount of data to produce them (mean = 0.96, SD = 0.01). This degree of density is important, as these methods are known to face challenges when computed on sparse networks (Fortunato and Hric 2016).

Figure 2 shows a stereotypical subdivision of the ROIs for a single subject (100307) using spectral partitioning. This method dissects a network by the positive and negative values derived from the Fiedler vector (the eigenvector related to the second-to-lowest eigenvalue of the network’s Laplacian matrix, see Methods), which we used to binarize the ROIs Since this is a data-driven approach, there is no a priori labeling for these communities. We thus defined the DMN community as the one containing the one DMN-specific area from the meta-analysis (7m), and ensured that the binary labeling of the vertices was consistent across subjects. In general, the resulting pattern shows large DMN coverage in PCC, with non-DMN vertices populating its anterior segment. mPFC, on the other hand, is mostly covered by DMN vertices in its ventral-anterior and dorsal-anterior areas, with a persistent non-DMN pocket between them. This non-DMN section extends posteriorly onto pregenual cinglate cortex. Note that vmPFC does not encompass area a25, which is often related to negative emotional valuation (Joyce and Barbas 2018). Interstingly, this topographical combination of dorsal and ventral mPFC, with extensive PCC coverage, matches previous meta-analytic patterns more closely attributed to social cognition than SV (Hiser and Koenigs 2018).

Figure 2. Brain partition for an example subject (100307). Fiedler vector values (top) are mapped onto the brain surface, dividing it into positive and negative communities. The bottom brain shows the binarized Fiedler vector, with yellow areas denoting the DMN community (as indicated by coverage of area 7m).

Before determining the degree of generazibility of this topological pattern across individuals, we examined the validity of these partitions by comparing them to patterns derived with a modularity algorithm (Clauset, Newman, and Moore 2004). Modularity seeks to find the pattern of communities that maximizes a quality function that compares the connectivity strength of any two nodes with a null probability model. Since modularity is agnostic to the desired number of clusters, it was capable of finding more than two communities in our datset. Thus, if both methods agreed in the number and topological organization of the current communities, it would add confidence to our results. We quantified the cross-method agreement with the adjusted rand index (ARI). ARI measures the proportion of pairs of nodes in a network that were either clustered together or separately in both partitions, over all possible affiliation changes (e.g. clustered together in partition 1, but separately in 2). This measure has the benefit of being agnostic to labeling schemes, as well as controlling for the expectation that the clusters were randomly produced (Hubert and Arabie 1985). This means that positive ARI values compare favorably against chance clustering. Figure 3 shows the distribution of ARI values comparing topologies produced by modularity and spectral partitioning across subjects.

Figure 3. Distribution of ARI indices when comparing modularity and spectral partitioning. The distribution shows overall high agreement. (Will probably remove this figure)

The distribution suggests that both clustering methods were in high agreement across subjects (mean = 0.85 sd = 0.11). The importance of this result lies in that while spectral partitioning forces a bipartition of the network, modularity does not. This adds evidence to the division of these areas into DMN and non-DMN during rest.

Next, we examined how similar the partitions were across subjects by computing the ARI between every pair of subjects, and found that there was considerable agreement overall (mean = 0.25, SD = 0.05). However, qualitative inspection of the brain community topologies shows good alignment for PCC, with a topologically similar mPFC pattern that shifts topographically across subjects. To quantify this heterogeneity in mPFC, we performed the same analysis for each region separately (Figure 4). Indeed, this analysis showed that the functional topology of PCC is well maintained across subjects (mean = 0.43, SD = 0.1), while mPFC (mean = 0.11, SD = 0.05) showed significantly lower levels of agreement (paired permutation, p = 0; Cohen’s D = 4.21).

Figure 4. Similarity matrix showing ARI values among all subjects for PCC (lower triangle) and mPFC (upper triangle) separately. PCC is well aligned across individuals, while mPFC shows clear topographical heterogeneity.

Partition Comparison: Day 1 vs Day 2

In order to determine whether the low agreement shown in mPFC across subjects was a function of intrinsic variability or a subject-specific pattern, we compared the agreement among communities estimated with data from each separate day. If the latter were true, within-subject ARI values should be consistently higher than across subjects. Since SP and modularity showed converging estimates, from here on we focus on results computed with spectral partitioning to ensure bipartitions. Communities were estimated on concatenated pairs of daily left-right and right-left phase encoded sessions (see Methods for details), and once again compared using ARI. Figure 5 shows pairwise comparisons among ten example subjects for PCC and mPFC separately.

Figure 5. Similarity matrices for 10 example participants showing pattern agreement across days for the whole partition (left), and for PCC and mPFC separately (right). These figures show high agreement within (but not between) subjects, with mPFC again showing significant variability.

Figure 6. mPFC displays higher within-subject agreement than between subjeccts, as evidenced by a high within/between subject mean ARI. PCC also shows some level of topographical idiosyncrasy (ratio > 1), but significantly lower than mPFC.

Once again we found greater alignment across subjects for PCC (mean = 0.39, SD = 0.1) than mPFC (mean = 0.07, SD = 0.04), but both areas showed higher levels of agreement within subjects (PCC: mean = 0.66, SD = 0.09; mPFC: mean = 0.29, SD = 0.1). However, a coarse comparison of mean ARIs neglects the degree to which subject-wise partitions were indeed idiosyncratic. Therefore, we divided each subject’s ARI between days by the ARI values with the days from the other subjects. This ratio denotes the relative specificity of their partition. That means, partitions that align well across subjects should yield ratios close to 1, while subject-consistent and unique partitions would produce high ratios. Figure 6 shows the ratios for PCC and mPFC for all subjects. As expected, we found that while both regions had good specificity (i.e. ratios > 1), the ratios for PCC were significantly lower than mPFC (paired permutation,p= 0; Cohen’s D = 2.9343989). These test/retest results suggest that the topographic varability seen in mPFC is reliably rooted in subject-specific organization, and not an inability to capture its topology.

Partition Comparison: Session

We extended the analysis of daily data by examining whether the organization of the DMN could be extracted using session data only. It is worth noting that the length of each session (approx. 14 mins long) falls below the stability threshold for fMRI-based modularity estimations (Gordon et al. 2017). Nonetheless, high ARI ratios with reduced data could signify the ability to localize a subject’s DMN organization with fMRI sessions that are easier to fit along with other study sessions of interest.

Figure 7. Similarity matrices showing pattern agreement across sessions for the whole partition (left), and for PCC and mPFC separately (right) (10 example subjects). The pattern from the daily analysis is preserved, with marked individuality for mPFC, but not PCC.

Computing SP on networks from individual sessions shows that we can still coarsely capture a subject’s unique organizational pattern (Figure 7), even though the overall levels of agreement decreased (PCC between subjects: mean = 0.31, SD = 0.11; PCC within subjects: mean = 0.51, SD = 0.14); mPFC between subjects: mean = 0.03, SD = 0.03; mPFC within subjects: mean = 0.14, SD = 0.08). We again computed each subject’s ARI ratio in order to quantify the specificity of these partitions, but this time using the mean ARI acros sessions (Figure 8).

Figure 8. Whole group ratio between within- and between-subject ARI for mPFC and PCC based on sessions. Consistent with the daily analysis, both regions show idiosyncrasy, with mPFC being significantly more subject-specific.

We can see again that the ratios for mPFC were higher than those of PCC (permutation,p= 0; Cohen’s D = 2.66), but both were mostly above 1. This further confirms that the intrinsic functional organization of mPFC is uniquely arranged per individual, and provides evidence that such patterns can be identified with relatively small amounts of data in spite of a reduction in the estimates.

Pattern stability over time

Figure 9. ARI between communities computed at each sliding epoch and the original partition for all subjects (left), as well as the distribution of ARI standard deviations per subject (right).

While the distinction of DMN and non-DMN is consistent across session and days, it is possible that the decrease in global ARI values going from full to session-specific data could be due to a dynamic shift in DMN organization. In order to address this possibility, we performed a sliding window analysis on each subject’s full time series (20 min windows shifting at 1 min), comparing the community arrangement derived from each window with the original pattern estimated using the whole time series (Figure 9). Note that the length of the windows is well equipped to find stable community features in brain networks (Gordon et al. 2017). This analysis revealed persistently high levels of agreement for most subjects, with overall low variability (overall mean of each subject’s average ARI = 0.65; SD = 0.05).

We additionally performed this analysis on mPFC and PCC separately to determine whether these results were region-specific (Figure 10). Even though the within-subject variance remained low for both regions, we found that the average agreement per subject was significantly higher for PCC (paired permutation,p= 0; Cohen’s D = 3.67). Nonetheless, the mean ARI across each subject’s time series was still relatively high for both PCC (mean = 0.77; SD = 0.06) and mPFC (mean = 0.51; SD = 0.08), thus suggesting that the each individual’s DMN distribution is cosiderably stable in time.

Figure 10. ARI between each window and the original pattern for PCC and mPFC separately (left), and the mean per-subject ARI for the whole group (middle). Regardless of the mean ARI, standard deviation remained low for both regions (right).

We complemented these analyses by asking whether the magnitude of values in the Fiedler vector was associated with the stability of nodes throughout a subject’s time series (PERSONAL NOTE: It might be better to justify this from the use of this vector in statistical mechanics rather than mere curiosity, i.e. extreme FV values are given to elements at the peaks of the wave. Find a simple citation for this.). Specifically, we examined whether Fiedler vector values could predict the proportion of DMN affiliations per node by performing a logistic mixed effects model, with subject acting as a random effect (Figure 11). The model yielded a positive significant relationship between these features (\(\beta\) = 757.38, SE = 1.68,p< 0.001), meaning that vertices with higher absolute Fiedler vector values were more persistent in their relationship with their corresponding community in time. These analysis show that the original pattern identified with the full dataset is temporally stable in time, and that there is potential value in thresholding the Fiedler vector as a means to capture idiosyncratic and reliable DMN vertices in relatively short fMRI sessions (i.e. less than 20 mins).

Figure 11. Relationship between the magnitude of Fiedler vector values and the proportion of affiliations of a vertex with DMN throughout 37 sliding windows. Grey lines denote each subject’s association, while the black line shows the fit from a mixed effects logistic regression.

Potential Noise Effects (this part is challenging)

Figure 12. Mean tSNR for mPFC and PCC across participants. This measure was derived from the minimally preprocessed fMRI data from the HCP, and thus shows the base level of reltive noise in each region.

Even though our findings point to a stable representation of each subject’s unique DMN organization, it is still possible that our results were affected by considerable levels of signal noise often found in mPFC (its ventral segment in particular). As a first step to account for this, we computed the temporal signal to noise ratio (tSNR) for each vertex in each subject’s minimally preprocessed fMRI data (i.e. mean of each vertex’s time points over their standard deviation). The tSNR of all vertices were averaged per region and subject (Figure 12), which showed that indeed mPFC (mean = 63.2; SD = 7.52) was overall noisier than PCC (mean = 88.93; SD = 9.13) across subjects (permutation, p = 0, Cohen’s D = 3.08). To more precisely examine this difference, we performed subject-wise comparisons of the tSNR values of all vertices between PCC and mPFC (total 100 permutation tests, 5000 iterations each). About 99 percent of the comparisons proved significant at alpha = 0.05 (FDR corrected; mean Cohen’s D = 0.83). In addition, we found that the levels of tSNR for each vertex showed significant negative correlations with flexibility, a graph-theoretic index that measures the proportion of times a node in a network changed affiliation across partitions (mean r = -0.23). These findings indicate that, for any given subject, vertices with low SNR resulted in more variable participation between the communities, and these vertices were more common in mPFC. This could explain why mPFC previously showed lower stability levels.

Figure 14. Differences in mean R-squared values for each region per participant. The figure shows that PCC vertices tend to be more coupled with the overall dynamics of DMN across subjects.

To further interrogate the potential effect of noise in our results, we computed the variance explained by activity from each DMN vertex on the mean time series of the DMN community. If the signal quality issues persisted after preprocessing the data, we expected that noisy vertices would have the least explicative power on the activity of the DMN community as a whole. This also implied lower \(R^2\) values for mPFC. Indeed, we found that for each subject, \(R^2\) values tended to be higher for PCC (100 percent of Wilcoxon ranksum tests across subjects had p values < 0.001) (Figure 14). These differences had moderately high effect sizes (mean Cohen’s D = 0.98). These results suggest that, overall, vertices in mPFC tended to be less in synchrony with the activity of the DMN community as a whole. Furthermore, the degree tSNR showed a moderate correlation with \(R^2\) values across participants (mean r = 0.35; SD = 0.11), meaning that noisier vertices were less in synchrony with the general activity of the DMN.

Given this trend, we asked whether there is a “gradient of noise” along the z-axis of mPFC, such that more ventral vertices would both be noisier and explain the least amount of DMN variance. Here we found that while \(R^2\) are somewhat correlated with dorsal-ventral placement across subjects (mean r = 0.3; SD = 0.12), there is no clear relationship with the z-axis across subjects (Figure 15). Specifically, low \(R^2\) are distributed evenly along this axis, denoting a degree of orthogonality between these features. Conversely, tSNR showed a strong linear relationship with the position of the vertex along this axis (mean r = 0.86; SD = 0.05). These results contextualize the degree of synchrony of mPFC vertices in terms of potential noise. While still an important factor, it is unlikely that the degree of within-subject pattern stability observed here could be due mostly to highly structured mPFC noise.

Figure 15. Vertex tSNR values across all subjects are correlated with their position along the z-axis in mPFC, with vmPFC showing the lowest tSNR (left). However, R-square values do not show this strong relationship (right), suggesting that the coupling of each vertex is not necessarily influenced by signal noise.

Correlation vs community detection in mPFC

Next, it is worth examining what we gain from using graph theoretic approaches over a simple seed-based functional connectivity analysis from PCC. For this purpose, we computed the mean time series of all the vertices contained in the DMN PCC ROI from Yeo et al. (2011) parcellation, and correlated it with the activity of every mPFC vertex from our ROIs. The similarity matrix below shows Day 1 vs Day 2 comparisons for mPFC partitions generated with seed-based correlations (top-left quadrant) vs community detection (bottom-right) across all participants, as well as how much did these methods agree with each other (not much, but enough, which is good). The important features here are that both approaches are very consistent, but community detection seems slightly more stable across days (quantified in the next section). Also, the cross-method agreement is fair, meaning that both approaches are showing similar topological features, but that these are more stable using graph theoretic approaches.

Figure 16. Correlation matrix comparing the inter-day spatial stability of maps derived from seed-based functional connectivity (7m-based) and the Fiedler vector. The top-left quadrant represents seed-based FC maps, and the bottom-right the Fiedler vector.

The previous plot shows the inter-methodological correspondence with an example subjects. The final plot shows the spatial correlation between days for each method across all participants. Community reflects the data from the lower right quadrant, correlation belongs to the upper left, and both is the cross-approach correspondence (top right). Community detection had a high and significantly higher RI across days than a simple correlation (permutation p-value < 0.001; Cohen’s D = 1.14). This is important because it shows that we can derive very stable topological features using community detection, even though mPFC has relatively low signal to noise ratio. Low intra-subject and high inter-subject variability of DMN partitions is in agreement with previous work [T. O. Laumann et al. (2015); Gordon et al. (2017); Braga and Buckner (2017); Kong et al. (2018)). Graph-based community detection might thus be a better way to detect stable functional topologies, since high stability throughout days is a consistent feature of DMN topologies.

Figure 17. Day 1 vs Day 2 oefficients of spatial correlations for every method, as well as between methods. Community detection through spectral partitioning provides the most stable inter-day estimates, even though both methods show good level of agreement.

Structural correspondence

Finally, we examined whether these topological features were defined by the structure of the brain by correlating each subject’s Fiedler vector with their structural MRI data. We found that the Fiedler vector is uncorrelated with myelin density (mean correlation = -0.01; SD = 0.06), cortical folding (mean correlation = -0.02; SD = 0.04), and surface thickness (mean correlation = 0.11; SD = 0.07). While these results suggest that each subject’s functional topology is mostly independent of brain structure, we will perform more detailed analyses of this in the future (i.e. are there any outstanding structural features in a subject that match their unique topology?).

Discussion

For clinical implications, see Hiser & Koenigs 2018. In short, being able to identify functional idiosyncrasies without much data can yield important insights in diagnosis (meh, we always say that..) of a variety of psychological diseases associated with mPFC. It would be more interesting to identify the repercussions to damage in mPFC per individual. They mention that subgenual vmPFC is a common target for deep brain stimulation (or TMS) as a depression treatment. The refinement of the technique is thus coupled to being able to identify personalized functional topologies. They list a number of results related to mPFC functional connectivity too. In the end, their suggestions for clinical implications of mPFC studies rely on the individual measure of functional expanse. (note that social-relevant areas in their meta analysis look similar to our patterns, and the emotion part to Barbas latest)

It will also be good to note that there is concurrent debate on the ground truth of community detection. In other fields, identifying a community structure is not necessarily associated with the true organization of the network. What happens is that you might find a tangential grouping that gives you information different from your original goal (Clauset, Newman, and Moore 2004)

Mention the lack of interhemispheric distinction, citing Mackey and Petrides in regards to L-R diffs.

It’s also worth mentioning that PCC is largely covered by DMN. Then why does this area overlap with SV? I think the Mueller paper mentions that PCC is an important hub of the brain, and perhaps subserves many cognitive functions (recent Wash U paper on this too).

And write a portion on the lack of subcortical reward-related regions (i.e. striatum). Say something like ‘incorporating the ventral striatum did notrefine these partitions through either method. This could be due to the lack of signal modulation of this area at rest, and should be examined through behavioral manipulations of value in the future. Instead, here we focus on capturing the intrinsic functional heterogeneity of mPFC as a way to pave the path to refined analyses aimed at disentangling DMN and SV’.

Talk about how there are other community detection algorithms, and that modularity and SP over and underfit, respectively. However, SP eigenmaps have been found to provide relevant functional information. In addition, both modularity and SP rely on connectivity strength, and thus is amenable to the correlative nature of functional connectivity.

Add a portion showing that summing up the binarized maps produces a map similar to the one produced by Margulies et al. and others.

It might also be worth mentioning that previous work has established the shift in network types during rest through hidden Markov models (Vidaurre et al., 2016?). While this suggests that we should see fluctuations in affiliation, and not stability, Vidaurre’s work establishes the dynamics of prevalent function among networks, and not their topologies. Our work focuses on the topographical arrangement of the DMN network, which can be active or not while maintaining its form. So, the two works are compatible. Note, I played a bit with AR and MA processes per subject to try to understand the autocorrelation in the sliding window. First, both acf/pacf and arima (order = n,0,n) R functions suggested that the data was better fit by AR. Moreover, ARs from 1 to 10 had neglegible differences, so I stuck with 1 results for simplicity. As expected the data is significantly autocorrelated, and when I checked plotted the residuals they were mostly within the 0.15 bound (normally distributed). Coefficients ranged from 0.6 to 1. So, there is clear rythmicity in the agreement between each window and the overall pattern. However, the intercepts from the ARn processes gave more or less the same estimates as the mean rand index over all the windows. This is mostly expected, as we saw that the standard deviation of the rand index was low across participants. I am currently interpreting this as a signal that there is clear periodicity that is otherwise harmless for our results: in the end, the agreement remains stable throughout the 4800 TRs.

In terms of noise, there are a handful of things to keep in mind. First is that yes, mPFC is noisier overall, but that did not preclude the relatively high levels of agreement across days and sessions per subject, even though ARI heavily penalizes discrepancies (less than linearly). Second, low \(R^2\) values were found all throughout the z-axis, meaning that variance explained is not strictly determined by the cardinal position of the vertex (and while the correlation was significant due to the sheer number of items, the coefficient was relatively low). Right now this might be ok as a way to bypass noise comments, but it would also be interesting to take advantage of the Fiedler vector here. Since it’s linearly related to \(R^2\) per subject, perhaps finding a way to threshold these values would provide cleaner and even more consistent maps.

Acknowledgements

Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

Supplemental Materials

Table 1. Significant brain regions per literature

Literature	Significant Brain Regions
Default Mode Network	V1, V2, V3, 4, 3b, RSC, POS2, PCV, 7m, POS1, 23d, v23ab, d23ab, 31pv, 5m, 5mv, 23c, 2, 3a, a24, d32, p32, 10r, 9m, 10v, a10p, 10pp, OFC, OP4, PFt, DVT, PGp, PFop, PF, PFm, PGi, PGs, 31a, s32, p10p, p24, V6, V4, A1, 8Av, 8Ad, 9-46d, 31pd, MBelt, LBelt
Subjective Value	RSC, d23ab, a24, d32, p32, 10r, 9m, 10v, OFC, MI, Pir, AAIC, 25, s32, p24, V2, PSL, 31pv, 8BM, 10d, 45, 47s, FOP4, PF, PGi, 31a, pOFC

References

Acikalin, M. Yavuz, Krzysztof J. Gorgolewski, and Russell A. Poldrack. 2017. “A coordinate-based meta-analysis of overlaps in regional specialization and functional connectivity across subjective value and default mode networks.” Frontiers in Neuroscience 11 (JAN): 1–11. doi:10.3389/fnins.2017.00001.

Bartra, Oscar, Joseph T. McGuire, and Joseph W. Kable. 2013. “The valuation system: A coordinate-based meta-analysis of BOLD fMRI experiments examining neural correlates of subjective value.” NeuroImage 76. Elsevier Inc.: 412–27. doi:10.1016/j.neuroimage.2013.02.063.

Bassett, Danielle S, Perry Zurn, and Joshua I Gold. 2018. “On the nature and use of models in network neuroscience.” Nature Reviews Neuroscience. Springer US, 1–13. doi:10.1038/s41583-018-0038-8.

Belkin, Mikhail, and Partha Niyogi. 2003. “Laplacian eigenmaps for dimensionalitu reduction and data representation” 1396: 1373–96.

Braga, Rodrigo M., and Randy L. Buckner. 2017. “Parallel Interdigitated Distributed Networks within the Individual Estimated by Intrinsic Functional Connectivity.” Neuron 95 (2). Elsevier Inc.: 457–471.e5. doi:10.1016/j.neuron.2017.06.038.

Burgess, G. C., S. Kandala, D. Nolan, T. O. Laumann, J. D. Power, B. Adeyemo, M.P. Harms, S.E. Petersen, and D. M. & Barch. 2016. “Evaluation of Denoising Strategies to Address Motion-Correlated Artifacts in Resting-State Functional Magnetic Resonance Imaging Data from the Human Connectome Project.” Brain Connectivity 6 (9): 669–80.

Clauset, Aaron, M. E. J. Newman, and Cristopher Moore. 2004. “Finding community structure in very large networks” 066111: 1–6. doi:10.1103/PhysRevE.70.066111.

Clithero, John A., and Antonio Rangel. 2014. “Informatic parcellation of the network involved in the computation of subjective value.” Social Cognitive and Affective Neuroscience 9 (9): 1289–1302. doi:10.1093/scan/nst106.

Fiedler, Miroslav. 1975. “A Property of Eigenvectors of Nonnegative Symmetric Matrices and its Application to Graph Theory.” Czechoslovak Mathematical Journal 25 (100): 619–33.

Fortunato, Santo, and Darko Hric. 2016. “Community detection in networks: A user guide.” Physics Reports 659. Elsevier B.V.: 1–44. doi:10.1016/j.physrep.2016.09.002.

Fox, M. D., A. Z. Snyder, J. L. Vincent, M. Corbetta, D. C. Van Essen, and M. E. Raichle. 2005. “From The Cover: The human brain is intrinsically organized into dynamic, anticorrelated functional networks.” Proceedings of the National Academy of Sciences 102 (27): 9673–8. doi:10.1073/pnas.0504136102.

Garcia, Javier O., Arian Ashourvan, Sarah Muldoon, Jean M. Vettel, and Danielle S. Bassett. 2017. “Applications of community detection techniques to brain graphs: Algorithmic considerations and implications for neural function.” bioRxiv, 209429. doi:10.1101/209429.

Ghasemian, Amir, Homa Hosseinmardi, and Aaron Clauset. 2018. “Evaluating Overfit and Underfit in Models of Network Community Structure,” 1–17. doi:arXiv:1802.10582v2.

Glasser, Matthew F., Timothy S. Coalson, Emma C. Robinson, Carl D. Hacker, John Harwell, Essa Yacoub, Kamil Ugurbil, et al. 2016. “A multi-modal parcellation of human cerebral cortex.” Nature 536 (7615). Nature Publishing Group: 171–78. doi:10.1038/nature18933.

Glasser, Matthew F., Stamatios N. Sotiropoulos, J. Anthony Wilson, Timothy S. Coalson, Bruce Fischl, Jesper L. Andersson, Junqian Xu, et al. 2013. “The minimal preprocessing pipelines for the Human Connectome Project.” NeuroImage 80. Elsevier Inc.: 105–24. doi:10.1016/j.neuroimage.2013.04.127.

Gordon, Evan M., Timothy O. Laumann, Adrian W. Gilmore, Dillan J. Newbold, Deanna J. Greene, Jeffrey J. Berg, Mario Ortega, et al. 2017. “Precision Functional Mapping of Individual Human Brains.” Neuron 95 (4). Elsevier Inc.: 791–807.e7. doi:10.1016/j.neuron.2017.07.011.

Gratton, Caterina, Timothy O. Laumann, Ashley N. Nielsen, Deanna J. Greene, Evan M. Gordon, Adrian W. Gilmore, Steven M. Nelson, et al. 2018. “Functional Brain Networks Are Dominated by Stable Group and Individual Factors, Not Cognitive or Daily Variation.” Neuron. Elsevier Inc., 439–52. doi:10.1016/j.neuron.2018.03.035.

Hiser, Jaryd, and Michael Koenigs. 2018. “The Multifaceted Role of the Ventromedial Prefrontal Cortex in Emotion, Decision Making, Social Cognition, and Psychopathology.” Biological Psychiatry 83 (8). Elsevier Inc: 638–47. doi:10.1016/j.biopsych.2017.10.030.

Hubert, Lawrence, and Phipps Arabie. 1985. “Comparing partitions.” Journal of Classification 2 (1): 193–218. doi:10.1007/BF01908075.

Joyce, Mary Kate P., and Helen Barbas. 2018. “Cortical connections position primate area 25 as a keystone for interoception, emotion, and memory.” The Journal of Neuroscience 38 (7): 2363–17. doi:10.1523/JNEUROSCI.2363-17.2017.

Kable, Joseph W., and Paul W. Glimcher. 2007. “The neural correlates of subjective value during intertemporal choice.” Nature Neuroscience 10 (12): 1625–33.

Kong, Ru, Jingwei Li, Csaba Orban, Mert Rory Sabuncu, Hesheng Liu, Alexander Schaefer, Nanbo Sun, et al. 2018. “Spatial Topography of Individual-Specific Cortical Networks Predicts Human Cognition, Personality and Emotion.” bioRxiv 36 (24): 213041. doi:10.1523/JNEUROSCI.4402-15.2016.

Kragel, Philip A., Michiko Kano, Lukas Van Oudenhove, Huynh Giao Ly, Patrick Dupont, Amandine Rubio, Chantal Delon-Martin, et al. 2018. “Generalizable representations of pain, cognitive control, and negative emotion in medial frontal cortex.” Nature Neuroscience. Springer US, 1. doi:10.1038/s41593-017-0051-7.

Laird, A. R., S. B. Eickhoff, K. Li, D. A. Robin, D. C. Glahn, and P. T. Fox. 2009. “Investigating the Functional Heterogeneity of the Default Mode Network Using Coordinate-Based Meta-Analytic Modeling.” Journal of Neuroscience 29 (46): 14496–14505. doi:10.1523/JNEUROSCI.4004-09.2009.

Laumann, Timothy O., Evan M. Gordon, Babatunde Adeyemo, Abraham Z. Snyder, Sung Jun Joo, Mei Yen Chen, Adrian W. Gilmore, et al. 2015. “Functional System and Areal Organization of a Highly Sampled Individual Human Brain.” Neuron 87 (3). Elsevier Inc.: 658–71. doi:10.1016/j.neuron.2015.06.037.

Levy, I., S. C. Lazzaro, R. B. Rutledge, and P. W. Glimcher. 2011. “Choice from Non-Choice: Predicting Consumer Preferences from Blood Oxygenation Level-Dependent Signals Obtained during Passive Viewing.” Journal of Neuroscience 31 (1): 118–25. doi:10.1523/JNEUROSCI.3214-10.2011.

Logothetis, Nikos K. 2008. “What we can do and what we cannot do with fMRI.” doi:10.1038/nature06976.

Mueller, Sophia, Danhong Wang, Michael D. Fox, B. T.Thomas Yeo, Jorge Sepulcre, Mert R. Sabuncu, Rebecca Shafee, Jie Lu, and Hesheng Liu. 2013. “Individual Variability in Functional Connectivity Architecture of the Human Brain.” Neuron 77 (3). Elsevier Inc.: 586–95. doi:10.1016/j.neuron.2012.12.028.

Tian, Ye, and Andrew Zalesky. 2018. “Characterizing the functional connectivity diversity of the insula cortex: Subregions, diversity curves and behavior.” NeuroImage. Elsevier Inc. doi:10.1016/j.neuroimage.2018.08.055.

Tobyne, Sean M., David E. Osher, Samantha W. Michalka, and David C. Somers. 2017. “Sensory-biased attention networks in human lateral frontal cortex revealed by intrinsic functional connectivity.” NeuroImage 162 (August). Elsevier Ltd: 362–72. doi:10.1016/j.neuroimage.2017.08.020.

Van Essen, D. C., K. Ugurbil, E. Auerbach, D. Barch, T. E.J. Behrens, R. Bucholz, A. Chang, et al. 2012. “The Human Connectome Project: A data acquisition perspective.” NeuroImage 62 (4): 2222–31. doi:10.1016/j.neuroimage.2012.02.018.

Vega, Alejandro de la, Luke J. Chang, Marie T. Banich, Tor D. Wager, and Tal Yarkoni. 2016. “Large-Scale Meta-Analysis of Human Medial Frontal Cortex Reveals Tripartite Functional Organization.” The Journal of Neuroscience 36 (24): 6553–62. doi:10.1523/JNEUROSCI.4402-15.2016.

Yeo, B.T. Thomas, Fenna M. Krienen, Jorge Sepulcre, Mert R. Sabuncu, D. Lashkari, Marisa Hollinshead, Joshua L. Roffman, et al. 2011. “The organization of the human cerebral cortex estimated by intrinsic functional connectivity.” Journal of Neurophysiology 106: 1125–65. doi:10.1152/jn.00338.2011.

Zilles, Karl, Nicola Palomero-Gallagher, and Katrin Amunts. 2013. “Development of cortical folding during evolution and ontogeny.” Trends in Neurosciences 36 (5). Elsevier Ltd: 275–84. doi:10.1016/j.tins.2013.01.006.